********************************************************************************
********************************************************************************
*************************** COST-BENEFIT ANALYSES ******************************

clear all
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"

*** prism-related variables
gen rsq_pre = max if treated==0
bysort Household : egen rsquare_pre = max(rsq_pre)
gen rsq_post = max if treated==1
bysort Household : egen rsquare_post = max(rsq_post)

gen hd_pre = HDDbest if treated==0
bysort Household : egen hdd_pre = max(hd_pre)
gen hd_post = HDDbest if treated==1
bysort Household : egen hdd_post = max(hd_post)

* prism sample restriction
gen prismrest = 0 if rsquare_pre!=. & rsquare_post!=. & hdd_pre!=. & hdd_post!=.
replace prismrest = 1 if hdd_pre>=43 & hdd_pre<73 & hdd_post>=43 & hdd_post<73 & ///
		rsquare_pre>=0.85 & rsquare_pre!=. & rsquare_post>=0.85 & rsquare_post!=.
		
gen aux = 1
bysort Household : gen homeobs_prism = sum(aux) if regsample==1 & treated==1 & prismrest==1
replace homeobs_prism = . if homeobs_prism!=1
drop aux
		
** merge with info about furnace replacements
merge m:1 Household using "$maindir/Data/furn_replace.dta"
drop if _merge==2
drop _merge


***** some extra restrictions on number of homes served by each contractor
gen aux=1
sort Household meterend
bysort Household : gen homeobs = sum(aux) if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs = . if homeobs!=1

sort Household meterend
bysort Household : gen homeobs_post = sum(aux) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs_post = . if homeobs_post!=1

** fixing builddate to avoid omitted variables
replace builddate = 11 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

** fixing sqfeet to avoid omitted variables
drop roundsqfeet
gen roundsqfeet = .
	forval i = 0(500)2500 {
		replace roundsqfeet = `i' if sqfeet>`i' & sqfeet<=`i'+500
	}
replace roundsqfeet = 3000 if sqfeet>3000 & sqfeet!=.

** fixing some spending variables to avoid omitted bins
drop binAirCon
gen binAirCon = .
replace binAirCon=0 if Real_tot_actAirCon==0
replace binAirCon=1 if Real_tot_actAirCon>0 & Real_tot_actAirCon!=.

replace binAttic = 2700 if Real_tot_actAttic>2400 & Real_tot_actAttic!=.

** demographic variables
gen simpocc = noccupants
replace simpocc = 5 if noccupants>=5 & noccupants!=.

gen binAge = .
forval i = 20(10)80 {
replace binAge = `i' if Age>`i' & Age<=`i'+10
}

gen binIncome = .
forval i = 5000(5000)35000 {
replace binIncome = `i' if Real_income>`i' & Real_income<=`i'+5000
}
replace binIncome = 1 if Real_income<=5000
replace binIncome = 40000 if Real_income>40000 & Real_income!=.


*** housing structure
gen binSqft = .
forval i = 600(300)2700 {
replace binSqft = `i' if sqfeet>`i' & sqfeet<=`i'+300
}
replace binSqft = 1 if sqfeet<=600
replace binSqft = 3000 if sqfeet>3000 & sqfeet!=.

drop BlowerPreBins
gen BlowerPreBins = .
forval i = 2000(1000)9000 {
replace BlowerPreBins = `i' if Blower_Pre>`i' & Blower_Pre<=`i'+1000
}
replace BlowerPreBins = 1 if Blower_Pre<=2000
replace BlowerPreBins = 10000 if Blower_Pre>10000 & Blower_Pre!=.

*** energy prices
gen binelec = .
local counter = 2
forval i = 10.5(0.5)12.5 {
replace binelec = `counter' if elec_prices>`i' & elec_prices<=`i'+0.5
local counter = `counter' + 1
}
replace binelec = 1 if elec_prices<10.5
replace binelec = `counter' if elec_prices>13 & elec_prices!=.

gen bingas = .
local counter = 2
forval i = 7(1)16 {
replace bingas = `counter' if gas_prices>`i' & gas_prices<=`i'+0.5
local counter = `counter' + 1
}
replace bingas = 1 if gas_prices<7
replace bingas = `counter' if gas_prices>17 & gas_prices!=.

* transforming decimals in percent
gen savings_gap_pct = savings_gap*100

forval i = 1/200 {
gen savings_gap_pct`i' = 100*savings_gap`i'
}

** counting number of homes served by contractor
sort arch_contID
by arch_contID : egen totjobs = total(aux) if ww_treat==1
by arch_contID : egen totaljobs = max(totjobs)
drop totjobs
replace totaljobs=. if arch_contID==.
* contractor ID 236 to be ommitted, being one of the biggest contractors in this sample

** main heat pre-treatment
gen MainHeatkBTU = MainHeatBTU/1000

gen binheatkbtu = .
forval i = 60(10)120 {
replace binheatkbtu = `i' if MainHeatkBTU>`i' & MainHeatkBTU<=`i'+10
}
replace binheatkbtu = 1 if MainHeatkBTU<=60
replace binheatkbtu = 130 if MainHeatkBTU>130 & MainHeatkBTU!=.

gen heatworks = 0 if MainHeatOperational=="False"
replace heatworks = 1 if MainHeatOperational=="True"


* fixing variables that have a zero indicator (for later interacting)
foreach x in nstories binAirCon white black hispanic otherrace haselderly hasminor heatworks {
replace `x' = `x'+1
}
replace roundAtticR = 1 if roundAtticR==0

gen heatnotwork = 1 if MainHeatOperational=="False"
replace heatnotwork = 0 if MainHeatOperational=="True"

replace binFurnace = 1 if binFurnace==0

gen failprism = 0 if prismrest==1
replace failprism = 1 if prismrest==0

* variables for projected savings
gen ExistingConsumption_mmbtu = ExistingConsumption/1000000
gen ProjectedConsumption_mmbtu = ProjectedConsumption/1000000


********** savings per bins of amount spent
***** mmbtu savings
set matsize 10000		
reg tau_ml ib1.binFurnace ib1.binAirSeal ib1.binAttic ib1.binWallIns ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu ib2.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		ib1.white ib1.black ib1.hispanic ib1.otherrace ib1.haselderly ib1.hasminor ///
		ib1500.binSqft ib1.nstories ib1.roundAtticR ib1.builddate ///
		ib2011.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		if regsample==1  & treated==1	
est sto savefull
sca obs = e(N)
sca r2 = e(r2)
di "R-squared for regression on savings = `=scalar(r2)'"
predict fullsavepredicts if regsample==1  & treated==1	

*** regression excluding household demographics - just to obtain an r-squared reported in the paper
qui: reg tau_ml ib1.binFurnace ib1.binAirSeal ib1.binAttic ib1.binWallIns ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu ib2.heatworks ///
		ib1500.binSqft ib1.nstories ib1.roundAtticR ib1.builddate ///
		ib2011.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		if regsample==1  & treated==1	
sca r2 = e(r2)
di "R-squared for regression on savings (excluding demographics) = `=scalar(r2)'"


*** bootstrapping the savings regression
forval i = 1/200{
dis "Starting iteration `i'."

qui : reg  tau`i' ib1.binFurnace ib1.binAirSeal ib1.binAttic ib1.binWallIns ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu ib2.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		ib1.white ib1.black ib1.hispanic ib1.otherrace ib1.haselderly ib1.hasminor ///
		ib1500.binSqft ib1.nstories ib1.roundAtticR ib1.builddate ///
		ib2011.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		[fweight = weight_boots`i'] if regsample==1  & treated==1 
		
predict savepredicts`i' if regsample==1 & treated==1 & tau`i'!=.

dis "Iteration `i' completed."
}


**** now create average savings per month
sort Household end_month
by Household end_month : egen use_post = mean(pred_full) if regsample==1  & treated==1	
by Household end_month : egen monthuse_counterpost = max(use_post)
drop use_post

by Household end_month : egen save_post = mean(fullsavepredicts) if regsample==1  & treated==1	
by Household end_month : egen monthsave_post = max(save_post)
drop save_post


*** bootstrapping the average savings per month
sort Household end_month
forval i = 1/200 {
by Household end_month : egen use_post = mean(pred_boots`i') if regsample==1  & treated==1	
by Household end_month : egen monthuse_counterpost`i' = max(use_post)
drop use_post

by Household end_month : egen save_post = mean(savepredicts`i') if regsample==1  & treated==1	
by Household end_month : egen monthsave_post`i' = max(save_post)
drop save_post

dis "Iteration `i' completed."
}


save "$maindir/Results/Model_Outputs/cba_data.dta", replace
